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ABSTRACT 


The effect of rolling motion of a wing on the magnitude of 
error introduced due to the wing vibration when measuring atmo- 
spheric turbulence with a wind probe mounted on the wing tip is 
investigated. The wing considered has characteristics similar to 
that of a B-57 Cambera aircraft, and Von Karman's cross spectrum 
function is used to estimate the cross-correlation of atmospheric 
turbulence. Although the error calculated is found to be less than 
that calculated when only elastic bendings and vertical motions of 
the wing are considered, it is still relatively large in the frequency's 
range close to the natural frequencies of the wing. Therefore it 
is concluded that accelerometers mounted on the wing tip are needed 
to correct for this error, or the atmospheric velocity data must 
be appropriately filtered. 
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CHAPTER I 


INTRODUCTION 

The aim of this study is to determine the relative error in turbu- 
lent velocity measurements caused by the vibrations of the wing of a 
B-57 Cambera-type airplane when atmospheric turbulence is measured by a 
probe mounted on the wing tip. NASA is planning a program to measure 
turbulence data using a B-57 airplane. Whether or not accelerometers 
should be mounted on the wing tip in order to obtain accurate measure- 
ments is relevant to this program. Reference [1] contains a similar 
analysis. It uses spectrum analysis to measure the wing tip velocity 
spectrum and the error introduced in the turbulence measurement due to 
the wing vibration for the same aircraft. In the aforementioned work, 
the motion of the airplane was restricted to vertical motion and 
vertical bending of the wing. The restriction of the airplane to 
vertical motion and vertical bending while ignoring rolling, may cause 
an overestimation of the magnitude of the elastic modes. If rolling is 
taken into account, the magnitudes of the elastic modes will be less due 
to a relaxation of the elastic modes. In other words, the coupling 
between the rolling motion and the other types of motion (as will be 
seen later, the coupling is between the rolling mode and the antisym- 
metric modes only) may cause the wing to vibrate less than the amount 
predicted when the airplane is restricted to vertical motion, and hence 


♦Numbers in brackets correspond to similarly numbered references in 
the List of References. 
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the relative error in the measured turbulence would be less than that 
calculated in [1], 

A secondary aim of this work is to more fully document the computer 
program that was used in [1] with the modifications for rolling motion 
included. The complete computer code with annotations is given in the 
appendix. 

Before tackling the direct problem, a brief discussion of spectrum 
analysis is given. Since Wiener first published his classical monograph 
[2] in 1933, spectrum analysis has been increasingly used in wide and 
diverse scientific areas. That monograph showed for the first time how 
to use the Fourier integral as a link between two otherwise distinct 
branches of mathematics--namely, statistics and analysis. The complex 
form of the Fourier integral theorem states that if f(t) is absolutely 
integrable on the whole t axis and if f(t) is piecewise smooth on every 
finite interval, then the following equality holds: 


f(t) 


oo 



— 00 


f(x)e 1 “ (x ’ t) dx 


The function: 


( 1 . 1 ) 


F(to) 


/2iT 


f(A)e 1wX dA 


( 1 . 2 ) 


is called the Fourier transform of f(t). 
f (t) , then we have: 


If Equation 1.1 holds for 


f (t) 


1 

/2t t 


f F( w )e“ lwt da ) 


(1.3) 
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i.e., f(t) is the (inverse) Fourier transform of F(w). The advantage of 
transforming a function from one domain to another (e.g., from time 
domain to frequency domain) is that the complicated mathematical opera- 
tions on the original function (such as convolution and differentiation) 
are reduced to simple algebraic operations on the transformed function. 

A brief historical account of Fourier analysis can be found in [3]. 
Reference [4] is a highly readable book on Fourier analysis. 
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CHAPTER II 


ASSUMPTIONS ABOUT TURBULENCE AND AIRPLANE MOTION 

To calculate the response of an airplane to atmospheric turbulence, 
several elements are needed. These are: 

1. The statistical description of the turbulent field (the 
input. 

2. The calculation of the aerodynamic forces associated with 
the turbulent field (the gust forces). 

3. The calculation of the transfer functions which relate 
the airplane motion (or other quantities of interest) 
to the gust forces. 

4. The combination of the turbulence description with trans- 
fer function to obtain the output. 

Steps 2 through 4 depend strongly on the method adopted for 
describing the turbulent field in Step 1. The analysis that follows 
utilizes the two-point correlation function of the velocity component of 
interest. The turbulent field is considered to be homogeneous, isotro- 
pic, and momentarily frozen (Taylor's hypothesis). 

Another important factor in describing the atmospheric turbulence 
is its randomness in the flight path of the airplane. In the case of 
what is called one-dimensional gust, the gust field is considered to be 
random in the direction of the flight but is assumed uniform in the 
spanwise direction as depicted in Figure la. This assumption is satis- 
factory as long as the ratio of the span of the airplane to the turbu- 
lence length scale is small (less than one-tenth). However, if this is 
large, i.e., the turbulence length scale is small, the gust should also 
be considered random in the spanwise direction as depicted in Figure lb. 
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(a) One-dimensional gust (b) Two-dimensional gust 

Figure 1 Illustration of a one- and two-dimensional gust. 

Thus in the case of two-dimensional gusts, the problem of calculat- 
ing the response of an airplane is essentially one of determining the 
response of a linear system, i.e., the airplane, to a multi-dimensional 
stationary random process which in this case is the atmospheric turbu- 
lence. The assumptions of linearity and stationarity made in Reference 
[1] are also made in this study. 

The problem of calculating the response of an airplane to gust 
loads in the general case is a formidable one. To simplify the work 
required, some assumptions concerning the motion of the airplane are 
needed. The airplane is restricted to vertical motion and to distortion 
in the first few free-bending modes of the wing. In addition to 
vertical motion and vertical bending, the rolling of the wing about its 
symmetry axis is taken into consideration. The restriction of the air- 
plane to vertical motion and elastic bending of the wing is tantamount 
to the assumption that the moment of inertia of the fuselage about the 
x-axis (see Figure 2) is infinite. The inclusion of the rolling motion 
of the wing and ignoring the rolling motion of the fuselage is 
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Figure 2 Coordinate system. 

equivalent to the assumption that the moment of inertia of the fuselage 
is negligible. In this case the central part of the wing is taken to 
have high stiffness so that it can roll but not bend appreciably. The 
actual wing vibration is somewhere between these two extreme cases, 
i.e., the actual solution is in between that calculated in Reference [1] 
and the one calculated in this work. The motions and bendings of the 
wing considered are those caused by only the vertical component of the 
gust acting on the wing of the airplane. 

According to the preceding simplifications, the pitching and yawing 
motions of the airplane are ignored. Also, the contribution to roll by 
the gust acting on the tail is not included in the analysis. 

The airplane considered in this work is one that has wing charac- 
teristics similar to that of a B-57 wing. The system of axes employed 
is shown in Figure 2 where the spanwise coordinate, y, is the indepen- 
dent variable and the vertical displacement, w, is the dependent vari- 
able. The angle of rotation of the wing about the x-axis is denoted by 
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0 . The wing is assumed to be a flat plate having stiffness similar to 
a B-57 wing, and the aerodynamic forces appearing later are those 
calculated from strip theory. The distributions of the mass and that of 
the semi-chord across the semi -span are shown in Figures 3 and 4, 
respectively. For calculating the atmospheric turbulence power spec- 
trum, the von Karman cross spectral function as defined in [5] is 
assumed. 
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Figure 3 Distribution of mass across the semi -span of the wing. 
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Figure 4 Distribution of semi -chord across the semi -span of the wing. 
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CHAPTER III 


SPECTRUM ANALYSIS 


3.1 Output Spectrum of a Single Input 

Consider a random process P(t) and an ensemble member p(t) of this 
process. The autocorrelation function of P(t) is defined by: 

R p (t,t+x) = E[P(t)P(t+x)] (3.1) 

where E denotes the expected value. In this case the process P(t) is 
stationary and the autocorrelation function becomes a function only of 
time difference x; that is, 

R p (t,t+x) = R p (x) (3.2) 

Note that in practical situations we are usually forced to work 
with only one ensemble member of a process and, consequently, derive 
mean value, correlation function, etc. from this member. For this we 
assume the ergodic theorem, which allows time averages to be equated to 
the corresponding statistical averages, i.e.. 


E[P(t) ] = lim Yf 

T-*x> 


p(t)dt 


(3.3) 


R (x) = lim 

P t-** 


1 

2T 


p(t)p(t+x)dt 


(3.4) 


An even more useful statistical characteristic of a random process 
is its power spectrum which is defined as the transform of its 
correlation function: 
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f 


(3.5) 


(to) 


Rp(T)e" iwT dT 

— oo 


If p(t) is the input to a linear time invariant system, the 
response of this system is expressed as a convolution integral of p(t) 
and an influence function h(t): 


a(t) 


co 

r 


—oo 


p(A)h(t-A)dA 


(3.6) 


For the analysis represented in the next chapter, p(t) represents the 
gust at one station along the wing and h(s,,t) is the impulse response 
at the right wing tip. The Fourier transform of h(£,t) will be denoted 
as H(&,w) and will be defined as the right wing tip velocity due to a 
gust located at a station y-| along the wing. To simplify notation, 
h ( £ , t ) and H(&,w) will be written as h(t) and H(w) with the under- 
standing that they represent the impulse response and the transfer func- 
tion at the right wing tip. Now, the autocorrelation function of the 
output can be shown [6] to be related to the correlation function of the 
input as follows:. 


R ct (t) 


oo 

j R p (T+A r X 2 )h(x l )h(x 2 )dX l dX 2 

— oo 


(3.7) 


Taking the Fourier transform of both sides of Equation 3.7 yields: 


i> a (w) 


R e” la)T dx 
a 


(3.8a) 


where 
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CO 


00 


(3.8b) 


4> a (w) = 


h(x 1 ) 


h(x 2 ) 


R p (T+x r x 2 )e 


-1 U)T 


dxdx^dx^ 


The change of variable x = t + x-j - X 2 produces: 

4> a («) = jf h(x n )e 1 j h(x 2 )e 2 j R p e' lu,X dX (3.9) 

— 00 —CO —CO 

The first two integrals are recognized as H*(cu) and H(u>), respectively, 
and the third is <|> (w) as given by Equation 3.5. The preceding equation 
can be written: 

<f> a (“) = 4>p(^) | H(tu) | 2 (3.10) 

This equation is the basic equation used in spectrum analysis. The 
procedure to calculate the output spectrum, in the case of single input, 
is now clear: The transfer function H(w) is calculated first and then 

the square of its magnitude is multiplied by the assumed input power 
spectrum which, for the purpose of the analysis presented herein, is the 
von Karman spectrum function. 


3.2 Output Spectrum of Multiple Inputs 

Consider a linear system having multiple random inputs. In the 
case treated here, the wing of an airplane subjected to two-dimensional 
turbulence may be considered such a system. The variation of gust 
intensity in the spanwise direction produces different influence func- 
tions hU,y,t) at different sections of the wing. The total response 
then is the superposition of the response due to disturbances at each 
individual section. Hence, the response may be given by the following 
equation: 
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o(t) = 


(3.11) 


Z oo 
f 

P(y^)h(y,t-x)dxdy 

-z 

where z is the half length of the wing and p(y,x) is the input at time x 
at a station y with width dy. Also, it is understood that the response 
is being computed at z\ hence, z has been omitted from the parenthesis. 
An alternative form of Equation 3.11 is: 

z °° 

o(t) = h(y,X)p(y ,t-x)dxdy (3.12) 

-Z -00 

The correlation function of the output is: 

Rj(t) = a (t)a(t+T ) 

Z Z 0° °° 



-Z -Z -oo 

X dx 1 dx 2 dy 1 dy 2 (3.13) 

where the overbar denotes the time average. Expressed in terms of the 
correlation function, R , of the input. Equation 3.14 becomes: 

Z Z oo CO 

■ (■ 

R^(t) - h (y^ , X ^ ) h (y 2 , X 2 ) Rp (y -j ,y 2 , t+X^ -X 2 ) dx i dx 2 dy ^ dy 2 (3.14) 

-Z -Z -0° 

Taking the Fourier transform and then using the change of variable x = 
TfXi -X 2 , Equation 3.14 becomes: 

z z 

f 

4> (w) = H*(y 1 ,a))H(y 2 ,uj)4) (y 1 ,y 2 ,w)dy 1 dy 2 (3.15) 

J J ’ 

-z -z 

where <J> p (y-] ,y 2 .w) is a cross spectrum function. From the definition of 

the cross spectrum, it can be seen [6] that the following relation holds: 
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(3.16) 


4> p (y-| ,y 2 >w) = ♦p(y 2 »yi»“J 

where "*" denotes the complex conjugate. Because of Equation 3.16 and 
because of the symmetry of the domain of integration and by changing the 
limits of integration, the integrand (with new limits of integration) in 
Equation 3.15 can be written as: 

<J> p (y-| ,y 2 »oi)H(y 1 ,u>)H*(y 2 ,u>) + <i) p (y 2 ,y-| ,io)H(y 2 ,u)H*(y 1 ,«) 

= 4> p (y-, ,y 2 ,u , )H(y 1 ,u)H*(y 2 ,aj) + 4>*(y-| ,y 2 ,w)[H*(y 2 ,cu)H(y 1 ,«)]* 

= 4> p (y-|y 2 »cjL>)H(y-j ,w)H*(y 2 ,a)) + [<)>p(ypy 2 ,co)H(y-| ,w)H*(y 2 ,u))]* 

= 2Re[<fip(y 1 ,y 2 ,( U )H(y 1 ,a))H*(y 2 ,a))] (3.17) 

where Re[ ] is the real part of the quantity [ ]. Therefore, Equation 
3.15 becomes: 

£ £ 

» 

<t> a (o>) = 2Re[<j> p (y 1 ,y 2 ,a>)H(y 1 ,a))H*(y 2 ,co)]dy 1 dy 2 (3.18) 

-a y 2 

By the assumption of isotropy of the atmospheric turbulence, the input 
spectrum ^(y-] >y 2 » tu ) is a function only of the separation distance n, 
i.e., <j>p(ypy 2 ,w) = <f>p(y-]-y 2 ><*))• By this assumption and by making the 
change of variables n = y-j - y 2 and y = y 2 and by interchanging the 
order of integration. Equation 3.18 becomes: 

2£ £-n 

<t> ff (aj) = 4>p(n,w)2Re[ H(y ,co)H*(y+n ,oj)dydn] (3.19) 

0 -£ 

This equation relates the output spectrum <j> (to) to the cross spec- 
trum <)>p(n,u>) of the inputs. In this study, <}>p(n,w) is the cross spec- 
trum of atmospheric turbulence, and H(y,u>) is defined as the velocity of 
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the right wing due to a unit sinusoidal gust of frequency w located at y 
along the wing. Consequently, a(t), as defined in Equation 3.11, 
becomes the root mean square velocity of the right wing tip due to the 
gusts acting on the entire wing. Hence, <f> (u>) is the spectrum of the 
velocity of the right wing tip caused by the vertical components of the 
gust. 
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CHAPTER IV 


SPECTRUM OF WING TIP VELOCITY 

In the Chapter I, two types of gusts were considered: one- and 

two-dimensional gusts. These two cases are illustrated in Figure 1, 
page 5. For one-dimensional gusts, the intensity of the gusts is the 
same along the span, and hence there is no tendency to roll. But for 
two-dimensional gusts, the intensities of the gusts at different points 
of the span are different. These variations in gust intensities along 
the span give rise to a net rolling moment, which then results in roll- 
ing motion. 

Whether to consider a gust to be one- or two-dimensional depends on 
the relative sizes of the wing and the scale of turbulence. To elabo- 
rate on this statement, let us consider the gust to be composed of a 
number of sinusoidal gust components. The components with long wave- 
length tend to have uniform and small intensity along the span, and 
hence they produce no moment. The components with very small wavelength 
produce intensities along the wing that tend to cancel each other and 
again produce no tendency to roll. Figure 5 illustrates the effect of 
the components with long wavelength on the wing, and Figure 6 illus- 
trates the effect of the components with short wavelength on the wing. 

In contrast to the components with very long or very short wave- 
lengths, the components with intermediate wavelength are found to be the 
major contributors to the rolling motion. These intermediate wave- 
lengths are of the order from one to ten times the span of the wing. 
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Figure 5 Effect of gust with long wavelength on the wing. 



Figure 6 Effect of gust with short wavelength on the wing. 
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The effect of a sinusoidal with wavelength equal to that of the span is 
illustrated in Figure 7. 

The airplane considered in this study is a B-57 Cambera-type air- 
plane which has a wing span of 2 i = 66 feet. Houbolt [7] recommends a 
length scale of atmospheric turbulence of L = 300 feet, which is about 
1.8 times the span of the wing; a ratio that falls within the range of 
the intermediate components mentioned earlier. Hence, a two-dimensional 
spectrum analysis is required and the rolling motion should be taken 
into account. 

The questions concerning the statistical description of the turbu- 
lence being settled (recall the assumptions given in the Chapter I and 
the discussion in the last few paragraphs) and following the outline 



Figure 7 Effect of gust with wavelength equal to the length of the 
wing. 
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given in the Chapter I there still remains the task of calculating the 
aerodynamic forces associated with the turbulent field, the determina- 
tion of the transfer function from the equation of motion, and finally 
the calculation of the output spectrum. 

Let us first consider the free vibration of a beam. Consider an 
element of length dy. This element undergoes a vertical motion and 
rotary motion about its center of mass and a shearing deformation. The 
free-body diagram and the geometry for the beam element are shown in 
Figure 8. The following quantities can be defined: 
w = deflection of the centerline of the beam. 

-ry = slope of the centerline. 


\p = slope due to bending, 
w 
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If there is no shear deformation, the perpendicular to the face of the 
cross section will coincide with the centerline of the beam. In case of 
shearing deformation, these lines will be distinct and the shear angle 
will be equal to t|> - dw/dy. 

From flexure theory [8], we have these two elastic equations for 
the beam: 


' \> 


dw _ V 
dy ~ k'AG 


and 


(4,1) 


_ _M_ 
dy " El 


(4.2) 


where V is the shear force, M is the bending moment, A is the cross- 
sectional area, k' is a factor depending on the shape of the cross 
section, G is the modulus of elasticity in shear, and El is the bending 
stiffness. The equations of equilibrium of the moments and of the 
forces are given by: 



(4.3) 


and 



(4.4) 


where J and m are the rotary inertia and mass of the beam per unit 
length, respectively. 

Substituting Equations 4.1 and 4.2 into Equations 4.3 and 4.4, 
respectively, we get: 


jL(ft 
dy( ti dy. 


+ k'AG 


dw 

(dy 


- JiJj = 0 


and 


(4.5) 


18 



^ i { k ' AG (|r - *)} ' 0 


(4.6) 


which are the coupled equations of free vibration for the beam in the 
general case. Theoretically, these equations can be solved simulta- 
neously for arbitrary variations in elastic and inertial properties of 
the beam. In practice, closed solutions are difficult to obtain, 
except under certain simplifying assumptions. For example, if El and A 
are constant, these two equations can be reduced, after eliminating \p, 
to a single equation: 


El 


3 2 w 
3 y 2 


+ m 


3 z w 

dtJ 


- J + 


EIM 


k 'AG 


a 2 w Jm 3 4 w n 

3y 2 3t 2 k'AG 3t^ u 


(4.7) 


The wing we are dealing with can be considered as a slender beam in 
which the cross-sectional dimensions are small in comparison with the 
length, and for which rotary inertial effects and shear deformations may 
be neglected [9], In this case, the equations of motion. Equations 4.5 
and 4.6, can be reduced to one equation as follows: Since the shear 

angle is zero (equivalently, the shear modulus G is infinite), ^ should 
be equal to dw/dy, also the term Jip in Equation 4.3 is neglected, then 
Equations 4.1 and 4.3 can be written as: 

* ■ ay H-8) 

and 


V = 


dM 

37 


(4.9) 


Substituting the value of V from Equation 4.9 into Equation 4.4 we get: 


dfM 

dy- 


= -mw 


(4.10) 


The bending moment is related to the curvature by the flexure equation: 
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(4.11) 



As indi dated by Timoshinko, et al. [8], the error in calculating 
the natural frequencies of a beam with small cross-sectional dimension 
in comparison to its length using Equation 4.12 instead of Equations 4.5 
and 4.6 is very small. Therefore, Equation 4.12 is adequate for calcu- 
lating the natural frequencies of the beam, and it will be used in the 
analysis that will follow. 

Consider a wing with varying mass and stiffness along its span. If 
a sinusoidal gust is acting on a spanwise element of length, a, centered 
about a point, y-| , then the following equation for the balance of forces 
holds: 

|At(eI 0] = -mw - F„ + FgS(y ,y, ) (4.13) 


where 



is the force due to beam stiffness and mw is the inertial force as 
defined in the preceding discussion. is an aerodynamic force due to 
the motion of the wing and is given in [10] as: 

= irpb 2 w + 2irpVbC(k)w (4.14) 

where p is the density of the flight medium, b is the semi -chord of the 
wing, V is the mean flight speed, and C(k) is the Theodorsen function. 
The Theodorsen function is a function of the reduced frequency 
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(k = cob /V) of the motion. In terms of the Bessel functions, C(k) is: 

(4.15) 


„ V J 1 + v o> + V 1 (V 1 - V - <<V 0 + J , J 0> 

fj, T f 0 P * (Y, - J 0 >2 


All the functions appearing in Equation 4.15 are functions of the 
reduced frequency k. F R is the force due to the vertical gust and is 
given [9] by: 

F g = 2TT P V 2 b fe) K( k ) (4.16) 

. 

where K(k) is the Kussner function. K(k) can be expressed in terms of 
Theodorsen and Bessel functions: 

(4.17) 

The function 6(y 2 ,y-|) which acts on F R in Equation 4.13 selects the 
portion of the wing which is subjected to gust and is zero everywhere 
except between y^ - A/2 and y-j + a/2 where its value is unity. 

Substituting the forces defined by Equations 4.14 and 4.16 into 
Equation 4.13 we obtain the following differential equation: 

3 2 


K = C(J Q - i^) + i J-j 


3y2 


El 


3 2 W 


3y T J 


-mw - irpb 2 w - 2irpVbC(k)w + 2npV 2 bK(k) ^ 5 (y»yi) (4.18) 


The boundary conditions for a free-body wing are: 

w"(s,,t) = w 1 " (£,t) = w"(-£,t) = w 1 " ( -£ ,t) = 0 (4.19) 

These boundary conditions express the fact that there is no shear nor 
moment at the wing tips. 

Equation 4.18 is a forced-vibration equation for the wing, the 
free- vibration equation is: 

[EI(w"(y,t)]" = -mw(y,t) (4.20) 
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Assuming that there is a solution for the free-vibration equation of the 
form: 

w(y,t) = *(y)T(t) (4.21) 

we obtain, after substituting in Equation 4.20, 


1 d2T = JL 

T dt 2 m* 


[El*"]" 


(4.22) 


This equation is valid only if both sides are equal to some constant to 2 . 
Thus, Equation 4.22 is equivalent to the two equations: 


+ w 2 T = 0 (4.23) 

[El*"]" = mco 2 * (4.24) 

The initial conditions for Equation 4.23 are: 

T(0) = f (0) = 0 (4.25) 


The boundary conditions of Equation 4.24 are: 

(A) = <p (-A) = *"(A) = *"(-£) = 0 (4.26) 


which corresponds to the shear and bending moment being zero at the 
ends, as is the case for free ends. 

Physically, <j> represents the shape of a natural mode and w is the 
vibration frequency corresponding to this mode. There is an infinite 
number of values of w which satisfies Equation 4.20, and to each one of 
these values there corresponds a particular *. Thus, the solution of 
Equation 4.20 can be expressed as: 


w(y,t) = l *. (y)T. (t) (4.27) 

i=l 1 1 

The natural modes, *. , are orthogonal [9] because of the choice of the 
boundary conditions. The orthogonality condition is given by: 


22 



for i = j 
for i f j 


(4.28) 


■ 

l 

[El 4> i " ]" <f>jdy = 0,2 

mi> i <f>jdy = 


M i“i 


-a -a 

Substituting the value of w(y,t) from Equation 4.27 into Equation 
4.18 and using the orthogonality property of the modes we get a system 
of linear differential equations in terms of the natural modes: 


NIo^t 
m m m 


-M T - 
m m 


Trpb| 


N 

I 

j=l 


A mj T j 


N 

27 rpVb R C(k) l 

R J-l 


B .T. 
mj J 


o u ^ y l^ 

+ 2irpV 2 b R K( k) — y — <J> m (y 1 )a(y 1 )A 


(4.29) 


where A 


mj 


and B m j are aerodynamic cross terms given by: 



a2 Vj dy 


(4.30) 



-SL 

i 




(4.31) 


-l 

The term a(y) is the semi-chord distribution defined as a(y) = b(y)/b R . 

M m is the generalized mass of the mode $ corresponding to natural 
frequency w m as defined by Equation 4.28. M-| and are the mass and 

moment of inertia of the wing, respectively, and they are the general- 
ized masses of the rigid-body vertical and rolling motions. These rigid- 
body modes have frequencies “i = w 2 = the rigid-body rolling 

<t >2 was not included in the system of linear differential equations 
represented by Equation 4.29. As mentioned earlier, the restriction of 
the rigid-body motion to vertical motion only may cause an overestima- 
tion of the elastic vibrations of the wing. Hence, in the calculations 
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that will follow, $2 will be included in the system of linear diffential 
equations so that the effect of rolling is accounted for in the solution 
of the system of equations. 

Let us assume that the input and the output are sinusoidal having 
the forms: 

y=ue iks and T m = T m e iks (4.32) 

where the variables k and s are defined as follows: 

yj. 

k = — and s = jp- (4.33) 


Substituting the values of y and into Equation 4.29 and dividing 

j |/r 

throughout by irpV 2 S e , we obtain a system of linear algebraic 
equations: 


y f2 2 E = k 2 y E + 

m mm nrm 


N 

k 2 I 

j=l 


A mj C j + 


2ikC(k) 


N 

l 

j=l 


B mj £ j 


2b R K(k) 

+ — s — ^m ( yi )a( yi )A 


where 


(4.34) 



v m 7rpb R S 
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- = jn V 

Sn b D - 
R v 

The solution 1" of the system of linear algebraic equations. 

Equation 4.34, represents the amplitude of the modal response of the 

deflection of the wing to a sinusoidal gust. If b, = T, e^ s is differ- 

m m 

entiated with respect to time, the amplitude of the wing tip velocity of 
the modal response is obtained: 


T m V . T m 

St~- ,k bT= — 

K u y 


(4.35) 


The transfer function, which is the velocity of the right wing tip 
velocity due to a sinusoidal gust located at y-j , is then: 


N T.(y, ,w) 

H(y 1 ,«) = i l - a — : z — (4.36) 

j=3 y 

The summation in the preceding equation is carried over the bending 
modes only since the navigation system in the airplane can subtract the 
rigid-body motion from the turbulence data. 

The numerical procedure for calculating the velocity spectrum of 
the wing tip is outlined in the following paragraphs (for full details 
see Reference [1]). The complete computer code with annotations is 
given in the appendix. 

The numerical procedure is essentially divided into three sub- 
routines. The free-vibration problem is solved in the first subroutine. 
The forced-vibration problem is solved in the second subroutine. 

Finally, the velocity spectrum of the right wing tip is calculated in 
the third subroutine. 
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The first subroutine solves the free-vibration equation: 


[EI<f>"]" = mco 2 <{> for -i < y < a (4.37) 

with boundary conditions 

<P U) = <j 3 ,M (-4) = <t>"U) = = 0 (4.38) 

The coefficient m is the mass per unit length and El is the bending 
stiffness of the wing. Both m and El are a function of y. In [1], 
three different distributions of El were considered. These are classi- 
fied as: (1) standard wing, (2) stiff wing, and (3) flexible wing. The 

standard wing has a stiffness given by: 

El = 9 x 10 8 for y < |11 j 

El = 9 x 10 7 for y > | 11 | (4.39) 

The above values were determined using a static analysis assuming a 
maximum loading on the B-57 wing and a load factor of 10 g. The funda- 
mental frequency of this wing was found to be about 4 Hz. 

The stiff wing has a bending stiffness given by: 

El = 3 x 10 9 for y < | 1 1 | 

El = 3 x 10 8 for y > |11 | (4.40) 

These values were determined by a trial and error method so that the 
fundamental mode would have a frequency close to 7 Hz as determined 
from an in-flight experiment with the B-57 at NASA/Langley Research 
Center. 

The flexible wing was considered for comparison purposes; it has 
a bending stiffness: 

El = 9 x 10 6 for y < |11 1 

El = 9 x 10 5 for y > | 1 1 ] (4.41) 
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The stiff wing is the most representative of a B-57 wing, and 
therefore it is used in the numerical calculations of this study. For 
comparison purposes, the standard wing is also considered. 

El and m are defined symmetrically in the domain [-£,£] and so the 
differential equation is either symmetric or antisymmetric. Because of 
this, the problem needs to be solved only for 0 < y < z. In this case, 
the boundary conditions must be defined at the midspan of the wing. The 
conditions at the origin for a symmetric function are: 

+ '(0) = (0) = 0 (4.42) 


and the conditions for an antisymmetric function are: 


4>(0) = *"(0) = 0 (4.43) 

An initial value of cu is assumed, and then the fourth-order differ- 
ential equation is changed to a system of four first-order equations. A 
shooting method is used to estimate the complete set of initial condi- 
tions, after which a Runge-Kutta scheme is used to determine the solu- 
tion of the initial value problem. The value of the solution of the 
initial value problem at y = z is compared with the boundary conditions 
at y = a for the two-point boundary value problem. The estimate of the 
complete set of initial conditions is improved. The process is repeated 
until the initial value problem has the correct value for the boundary 
condition at y = a. The system of four first-order equations: 


_d_ 

dy 


*1 

/• 

<f >2 


4,3 


+4 

V 


0 

0 

m<u 2 

El 


1 

0 

0 

0 


0 

1 

0 

-2EI 

El 


0 

0 

1 

-El" 


El 


*1 

4>2 

$3 

<f>4 


(4.44) 
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of its solution. Each solution: 


has a vector base, , 

V 

*2 


(4.45) 


of Equation 4.43 is a linear combination of the base solution, i.e., 
4 

* = l ^i c i 

i=l 1 1 


For to be a complete base, there should be four linearly inde- 

pendent intial conditions. A simple choice of linearly independent 
initial conditions is: 


^1 (0) 


T 


'o' 


'o' 


o' 

0 

; ^o(o) = 

1 

; 'K(o) = 

0 

; 'Mo) = 

0 

0 

L 

0 

j 

1 

*+ 

0 

0 

v y 


0 

V ✓ 


lO. 


J. 


(4.46) 


The boundary conditions are used to determine the constants {C..} 
In the case of the even mode, we have the system of equations: 


^1 2 (°) ^ 22 ^ ^ 32 ^ ^ 42 ^°) 


fc 

L 1 


o' 

^14(0) ^24^) ^34(0) ^(O) 


C 2 


0 

^1 3 ( ^) ^ 23 ^ ^33 ( ^) ^43^) 


C 3 


0 

4 ( ^ ) ^24^^ ^34 ( ^ ) ^44 ( ) 




,0. 


4 

1=1 • 


(4.47) 


To have a nontrivial solution of the above equation, the determinant of 
the coefficient matrix must vanish. This determinant: 
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0 


1 


0 


0 


DU) = 


0 

^i 3 ( £ ) 

^ 14 (£) 


0 

ii> Z 3 (z) 
^24 ( ^ 


1 

^ 33 (£) 

* 34 (£) 


0 

^43 ( £ ) 
^44 ( ^ ) 


^U) 


^ 43 ( £ ) 

^44 ^ ^ ) 


(4.48) 


is a function of w, the unspecified parameter of the differential equa- 
tion. This characteristic determinant vanishes when the correct values 
of the natural frequencies, , of the wing are found. 

Now, we summarize the preceding numerical procedure. First, the 
eigenvalue is estimated; then the fundamental solutions are determined 
by a Runge-Kutta/Fehlberg order seven scheme. The value of the charac- 
teristic determinant is calculated from the fundamental solutions. A 
search routine checks if the natural frequency is bracketed between the 
current estimate and the previous estimate. In this case, the program 
is directed to a bisection routine to improve brackets or continues, 
taking another step along the frequency domain and using this as its 
next estimate of the frequency. After the frequency has been deter- 
mined, the natural mode is normalized by a unit displacement at the 
right wing tip. The new mode is integrated with previous modes deter- 
mined to calculate the aerodynamic cross terms. The program then steps 
along the frequency domain for its first and second estimates of the 
next frequency. 

The second subroutine of the program solves the system of linear 
algebraic equations. Equation 4.34. This system is solved by Gaussian 
elimination for 38 gust locations evenly spaced along the wing. 

The third subroutine calculates the output power spectrum. It 
integrates Equation 3.19 numerically using the trapezoidal rule: 
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N N-l 

<f> (n,w) = <f> (0,co) l H.H.* + l 4> (jA,iu)2Re 
a P j=l J J j=l P 


N-j 

l W 

i=l 1 3 


(4.49) 


where N is the number of gust stations and a is the gust station width. 
The spectrum program determines the wing tip velocity spectrum; there- 
fore, H.j must represent the velocity at the right wing tip due to a 
sinusoidal gust at station i. The frequency response function is 
defined as: 

NM 

H 1 ■ % j=l Tj (4-50) 

where NM is the number of elastic modes considered. The summation 
includes the response of only the elastic modes because, as mentioned 
earlier, the navigation system located at the airplane's center of 
gravity allows subtraction of the rigid body motions from the turbulence 
data taken at the wing tip. 
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CHAPTER V 


DISCUSSION AND CONCLUSIONS 

In the spectrum analysis section, a mathematical treatment of 
stationary processes was shown. In order to apply these techniques to 
problems dealing with atmospheric turbulence, the turbulence was assumed 
to be homogeneous. In general, most types of turbulence tend to be 
homogeneous, except those at low altitudes which may be influenced by 
the configuration of the ground and those in thunderstorms. 

The assumption of isotropy simplifies the expressions of two- 
dimensional correlation functions because in this case a correlation 
function can be expressed as one-dimensional correlation function of the 
separation distance. For sufficiently short wavelengths, all turbulence 
is isotropic [10], but for long wavelengths it is isotropic only if it 
is homogeneous. 

Taylor's hypothesis permits time correlation functions to be 
expressed as space correlation functions. This is so because the gust 
intensities remain the same until the airplane traverses this region of 
turbulence. The speed of the airplane plays an important role in 
validating this hypothesis. At very low speeds, Taylor's hypothesis 
becomes less valid and the results calculated may be less accurate. 

A quantity of great importance is the turbulence length scale L. 
Roughly speaking, L is a measure of the minimum separation distance for 
which there is no correlation between two velocity components of the 
gust. Many factors, such as weather conditions, altitude, and ground 
configurations (for low altitudes) play a role in determining L. The 
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estimates for L range between a few tens of feet to thousands of feet. 

In this study, two turbulence length scales of 132 feet and 2,112 feet 
were considered. Figures 9 and 10 show the turbulence spectrum and the 
wing tip velocity spectrum of a standard wing and a stiff wing for these 
two turbulence scales. 

The relative errors in the atmospheric spectrum due to wing tip 
vibration are shown in Figures 11 through 14. These were calculated 
using the following reasoning: The wind velocity measured at the wing 

tip is an apparent velocity since in reality it is the sum of the wing 
tip velocity and the true velocity of the wind. This can be expressed 
by the equation: 


U = U + U 
m p r 


(5.1) 


where U is the measured velocity, li is the true velocity of the turbu- 
m p 

lence, and U r is the wing tip velocity. From Equation 5.1 the correla- 
tion function of the measured velocity can be expressed by the following 
equation: 

T 


R m (x) = | U m (t)U m (t+T)dt 

T-*°° -T 


lim 

T-*» 


[U o (t)U a (t+T) + U p (t)U p (t+x) + U a (t)U p (t+x) 


+ U p (t)U a (t+x)]dt 

= R o (x) + R p (x) + R op (x) + R p0 (x) (5.2) 

Taking the Fourier transform of the above equation we get: 

♦m^) = i> a ( a3 ) + 'i’pU) + ^pU) + < f ) p a ( t °) 
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(5.3) 



Reduced Frequency (wb R /V) 

Figure 9 Spectra of wing tip velocity for standard and stiff wings 
in atmospheric turbulence having spectrum illustrated 
(length scale = 132 feet). 





Normalized Power Spectrum U and 4, ) 


Reduced Frequency (wb^/V) 


Figure 10 Spectra of wing tip velocity for standard and stiff wings 
in atmospheric turbulence having spectrum illustrated 
(length scale = 2,112 feet). 
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Reduced Frequency (wb^/V) 

Figure 11 Relative error in measuring turbulence with length 
scale 132 ft for standard wing. 




Relative Error [(|<j> 



Reduced Frequency (ub R /V) 


Figure 12 Relative error in measuring turbulence with length scale 
2,112 ft for standard wing. 


Relative Error [(|<j> m - 4> p l )/4> p 3 
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Reduced Frequency (cob R /V) 


Figure 13 Relative error in measuring turbulence with length scale 
132 ft for stiff wing. 
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Reduced Frequency (tub R /V) 

Figure 14 Relative error in measuring turbulence with length scale 
2,112 ft for stiff wing. 






(5.4) 


Using the relation <j> ap = <j>* a . Equation 5.3 becomes: 

<j> = <ji + 4> + 2Re<f> 
m a p T ap 

where cf> p is the turbulence spectrum at the wing tip and <j> 
spectrum of the wing tip velocity. The equation for 
as: 


is the power 
can be written 



2i 

* 

< f , D (n)H(£-n)dn 
J r 
0 


(5.5) 


The relative error due to the wing tip velocity can be written as: 



Figure 11 shows the relative error as a function of frequency for a 
standard wing and turbulence scale of 132 feet for two cases: (1) the 

rigid-body motion is restricted to vertical motion and (2) the rigid- 
body rolling is taken into account. This wing shows a large relative 
error close to its fundamental frequency. This error reaches a maximum 
of about 115 percent when the wing is restricted to vertical motion. 

This maximum is reduced to about 94 percent when the rolling motion is 
taken into account. The true error is believed to lie within these 
limits (see Chapter II, page 5). The range in which this maximum occurs 
falls within the frequency range which contains significant turbulence 
energy. For a turbulence length scale equal to 2,112 feet, the maximum 
relative error (see Figure 12) for the standard wing is reduced to about 
50 percent when rolling is not considered and to about 38 percent when 
rolling is added. But, this time the maximum relative error falls 
within a range which is out of the frequency range that contributes sig- 



nificant turbulence energy. The stiff wing shows a smaller relative 
error than the standard wing and the range of the maximum error is 
slightly shifted to higher frequency. The addition of rolling motion 
reduces this maximum from about 50 to 36 percent. For smaller turbu- 
lence scales, this maximum error is very close to the range of frequen- 
cies which contribute significant turbulence energy. But, for larger 
scales this error is outside the range of frequencies of interest. 

From the preceding discussion it can be seen that a stiff wing is 
the best wing to measure turbulence and but even then in order to measure 
accurately the whole range of the atmospheric spectrum depending on the 
length scale, accelerometers mounted on the wing tips are required. 
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APPENDIX 


DOCUMENTATION OF THE COMPUTER PROGRAM VIB.FOR 

The purpose of this appendix is to document the computer program 
VIB.FOR used in [1] and, with modifications, in the text. Flowcharts 
illustrating the complete computer code are given in [1]. An explana- 
tion of the subroutines and parameters appearing in the program will 
follow. 

VIB.FOR is essentially divided into three subprograms: (1) the 

free-vibration subroutine, D1 , (2) the forced-vibration subroutine, D2, 
and (3) the spectrum-analysis subroutine, D3. Each of D1 , D2, and D3 
has a number of subroutines and functions arranged in the following order 

I. SUBROUTINE D1 

1. Subroutine SZERO(W,H,N,NC,NN) 

2. Subroutine RUN(W,F,EL,EU,NC,NN) 

3. Subroutine BISEC(X1 ,F1 ,X2,F2,NC,NN) 

4. Subroutine FUNEV(K,X,Y,F) 

5. Function SEI(X) 

6. Subroutine SECANT(X1 ,F1 ,X2,F2,NC,NN) 

7. Subroutine DMODE(YY,NN) 

8. Function RM(X) 

9. Function SIMP(Y,H,N) 

10. Subroutine C0F(N,YY,W2,NN,I1 ) 

11. Subroutine STRK(NN) 

12. Subroutine RK7(NS,NN,EL,EU) 
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II. SUBROUTINE D2 


1. Subroutine D0(NN,N2) 

2. Subroutine C0EF(RK,X,DEL,N2M,NN) 

3. Function RPH(I,Y,N) 

4. Function RA(Y) 

5. Function RY1 (X) 

6. Complex Function CC(RK) 

7. Complex Function RKK(RK) 

8. Function RJ1 (X) 

9. Function RYO(X) 

10. Subroutine GAUSS (N,N2) 

11. Subroutine BACKS (N1 ,N2) 

III. SUBROUTINE D3 

1. Subroutine SPEC(RK,RR,TL,N2 s IC,RP) 

2. Function TSPEC(SS,RNU,U,TL) 

3. Subroutine C0EF1 

4. Function BSL1 (Z) 

5. Function P0LY(A,N,Z) 

6. Function BSL2 (Z ) 

7. Subroutine C0EF2 

A . 1 Subroutine D1 

In this subroutine the first four frequencies and elastic modes are 
calculated. Frequencies are calculated by starting with an initial 
guess after which the differential equation representing the free 
vibration is solved. From this solution the value of the characteristic 
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determinant is calculated. The first frequency is incremented and the 
differential equation is solved for this new value of the frequency. 

Then the value of the characteristic determinant is calculated and is 
compared with the previously calculated value. If these two values 
have the same sign, the frequency is incremented again and the differ- 
ential equation is solved and the characteristic determinant is calcu- 
lated. This process is done repeatedly until the value of the 
characteristic determinant changes sign, after which a bisection method 
and a secant method are used to approximate the frequency that makes 
the characteristic determinant vanish. The mode shape corresponding 
to this frequency is calculated, and then the generalized mass and the 
aerodynamic cross terms of all previously determined modes are calcu- 
lated. The same procedure is used until the first four natural frequen- 
cies and the corresponding elastic modes are determined together with 
all generalized masses and aerodynamic cross terms. 

Parameters . NC is the number of the mode to be calculated (the 
first value of NC is two because model number one is considered to be 
the rigid body translation mode which has zero frequency and a constant 
normalized value equal to unity. Modes number two, three, four, and 
five are the first four elastic modes. Mode number six is taken in 
this program to be the rigid body rolling mode). H is the step size 
in the frequency domain. W is the initial guess of the natural 
frequency. N = 4 is the number of elastic modes to be calculated. 

NN = 151 is the number of nodes on a semi -span at which the values of 
the modes are calculated. 
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COMMON/EDAT/EEE 
CALL D1 
CALL D 2 
CALL D3 
STOP 
END 

SUBROUTINE 111 

C*** THIS PROGRAM IS USED FOR DETERMING UING'S NATURAL 
FREQUENCY AND MODE. 


IT USES A BISECTION LIKE METHOD TO IMPROVE THE GUESS 
VANISH Ffi;E:QlJ£:NCY AND m AKE s THE CHARACTERISTIC DETERMINANT 

RUNGE-KUTTA FEHLBERG 7X8 IS USED TO DETERMINE SOLUTION 
TO DIFFERENTIAL EON IT HAS VAIRABLE SIZING BETWEEN 
FIXED NODES. 

IS USED T0 INTEGRATE FOR the areodynamic 
CROSS TERMS AND GENERALIZED MASS. 


THE PROGRAM MAYBE ADOPTED BY CHANGING FUN...SEI 
(THE FUNCTION DESCRIBING THE DISTRIBUTION OF THE PRODUCT 
YOUNG'S ELASTICITY AND SECTIONAL MASS MOMENT. IF SEI IS TO 
HAVE DERIVATIVES UP TO THE SECOND ORDER FUN . . . FUNEV NEEDS TO BE 
INCLUDE THEM IN THE DIFFERENTIAL EON. THE WING 
PLAN ARRAY. RA. UILL NEED TO BE CHANGED. 

C***NC=THE NUMBER OF MDDE TO START CALCULATING 
C***U=THE FRIST GUESS OF THE NATURAL FREQUENCY OF THE NC MODE 
C***H=STEP SIZE FOR SEARCH ROUTINE 
C***N=TOTAL NUHBER OF MODES TO CALCULATE 
C***NN=NUMBER OF NODES TO CALCULATE MODE ON 
REALMS U,H 
NC =2 
H=1 .5D0 
W=0 . IDO 
N = 4 

NN=I51 

CALL STRK(NN) 

CALL SZERO(W.H.N.NC.NN) 

RETURN 

END 


A. 2 Subroutine SZERO(W,H,N,NC ,NN) 

Description . This subroutine starts with an initial guess of the 
natural frequency and calls subroutine FUN which fixes the initial 
values of the differential equation. Then FUN calls subroutine RK7 
which solves the differential equation for the given frequency using 
Runge-Kutta-Fehlberg method of order seven. After obtaining the solu- 
tion of the differential equation, FUN calculates the value of the 
characteristic determinant and returns it to SZERO. Then, SZERO steps 
along the frequency domain and repeats the above procedure to calculate 
a second value of the characteristic determinant. SZERO checks for a 
change in sign of the calculated values of the determinant. If there is 
no change of sign, it steps again along the frequency domain and calcu- 
lates a new value of the characteristic determinant, and then it checks 
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for a change in sign of the last calculated values. If there is a 
change in sign, SZERO calls subroutine BISECTION to approximate the 
frequency that makes the characteristic determinant vanish. Then, SZERO 
calls subroutine SECANT to improve the approximation found by BISECTION, 
after which it calls subroutine DMODE to calculate the normalized values 
of the mode at 151 points of the semi-span. Finally, SZERO calls sub- 
routine COF to calculate the generalized mass of the mode and the aero- 
dynamic cross terms of the mode with the other modes. The frequency 
which was calculated last is used as an initial guess for the next 
frequency and the whole process is repeated until the first four fre- 
quencies are calculated. 

Parameters . I is a counter which gives the number of calculated 
modes. EL and EU are error bounds for Runge-Kutta method. W1 and W2 
are two consecutive values in the frequency domain. FI and F2 are the 
corresponding values of the characteristic determinant. 


c 

c 

c 

c 

c 

c 

c***« 
ctttt 
cttt* 
c***» 
ctttt 
c***» 
cttt* 
cttt t 
c***» 


103 

101 


102 


SUBROUTINE SZERO < W > H , N , NC . NN ) 

THIS SUBROUTINE STEPS ALONG THE FREQUENCY DOMAIN UNTILL 
IT BRACKETS THE NATURAL FREQUENCY AND THEN CALLS A BISECTION 
ROUTINE AND A SECANT ROUTINE TO IMPROVE THE UIDTH OF THE BRACKETS. 
THE SUBROUTINE FINALLY CALLS SUB... DMODE TO DETERMINE THE 
MODE r AND THEN CALLS COEF TO SET UP THE INTERGR AT I ONS OF MODES 
AND THEIR PRODUCTS. 

EL i EU=ERROR BOUNDS FOR RUNGE KUTTA 

VJ1 1 W2=BR ACKETS FOR FREQUENCY AND STEFS FOR THE SEARCH ROUTINE 
FI 1 F2=VALUES OF THE CHARACTERISTIC DETERMINTE FOR U1 t U2 
N=NUMBER OF EIGENVALUES TO BE SEARCHED 
U=GUESS FOR EIGEN VALUE 
H=STEP SIZE FOR SEARCH 

NC=NUMBER OF THE FRIST MODE TO CALCULATE 
NN=NUMBER OF NODES TO CALCULATE MODES ON 
•I =NUMBER OF MODES CALCULATED DURING PROCESS 
REAL*8 U.HiUl >W2>F1 > F2 > EL > EU > Y » X 
C0MM0N/FACT/Y<8f 151>,X<151),GM(6> 

DIMENSION Y Y ( 1 51 ) » XX (151) 

1=0 
U1 =w 

EL= . OOOIDO 
EU=.001D0 

CALL FUNIWl ,F1 . EL » EU . NC , NN ) 

U2=W1+H 

CALL FUN<W2,F2,EL»EU.NC.NN) 

IF<F1*F2.LT.O.O) GO TO 102 

W1 = U2 

F1 = F2 

GO TO 101 

CONTINUE 

CALL BISEC<W1»F1»W2»F2,NC»NN) 

CALL SECANTIUl rFl.U2.F2i NC,NN> 

CALL DMODE < YY » NN ) 

1 = 1 + 1 

US=U2 

CALL COF(NC,YY,USrNN,I) 

NC=NC+ 1 
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— IF ( I • GE • H ) GO TO 104 

U1=U2+H 
GO TO 103 

104 DO 8710 IjliNN 

071 a = _ 

CALL COF<A» YYf 0* »NN»5) 
RETURN 
END 


A . 3 Su broutine FUN(W,F ,EL,EU,NC,NN) 

Description . This subroutine estimates the initial conditions for 
the even and odd modes and then calls RK7 to solve the differential 
equation. Finally, it determines the characteristic determinant and 
returns it to SZERO. 

Parameters . EL and EU are error bounds for Runge-Kutta method. NC 
is the number of mode being determined. W is the frequency. F is the 
value of the characteristic determinant. Y is the matrix of the solu- 


tion of the system of differential equations solved in RK7. 


SUBROUTINE FUN ( U . F f EL . EU . NC , NN ) 

C THIS SUBROUTINE FIXES THE INITIAL VALUE FOR THE SOLUTION 

C AND THEN CALLS SUB . . RK7 < RUNGE KUTTA ROUTINE) TO DETERMINE SOLUTION 

C AFTER WHICH THIS SUB CALCULATES THE VALUE OF THE CHARACTERISTIC 

C DETERMINANT. 

C*** EL l EU =ERROR BOUNDS FOR RUNGE KUTTA 
C*** NC = UH I CH MODE WORKING ON DETERMING 

C*** Y=FUND AMENT AL SOLUTIONS TO DIFFERENTIAL EQN OUTPUT FROM RK7 
C»** F= VALUE OF CHARACTERISTIC DETERMINTE 
REAL*8 RKfWfFfEL.EUfYfX 
C0MM0N/FACT/Y<8fl51)fX(151) 

COMMON/EGNVL/RK 
RK = U 

C INITIAL CONDITION TEST FOR EVEN OR ODD MODE 
DO 100 1=1.8 

100 Y < I • 1 ) =0 • DO 

IF (NC.EQ.2.0R.NC .EQ. 4 ) GO TO 101 
C INITIAL CONDITION FOR ODD MODE 
Y<2» 1 >=1 . DO 
Y(8f 1 >=1 .DO 
GO TO 500 

C INITIAL CONDITION FOR EVEN MODE 

101 Y(1f1)=1.D0 
Y<7. 1 >=1 .DO 

500 CALL RK7(8fNN.ELfEU> 

F=Y<3fNN)*Y(8fNN)-Y<7fNN)*Y<4fNN> 

RETURN 

END 


A. 4 Subroutine B I SEC (XI ,F1 ,X2,F2,NC ,NN) 


Description . This subroutine approximates the frequency after 


finding an upper and a lower bound. 
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Parameters . XI is the lower bound of the frequency and X2 is the 


upper bound. FI and F2 are the values of the characteristic determinant. 


NC is the mode number of the frequency being determined. 


SUBROUTINE BISEC<X1fF1fX2fF2fNCf NN ) 
REAL*8 X1fX2fF1fF2fFHfXMfELfEU 
CCB= .01 
EU= . 001D0 
EL= • OOOIDO 

102 XM=(X1+X2)/2.D0 

CALL FUNIXMfFMfELfEUfNCfNN) 
IF<F1*FM.LE.O.BO) GO TO 100 


GO TO 101 

100 X2=Xh 
F2=FM 

101 RE = BABS < X1-X2 ) /XI 
IF<RE.GT . CCB > GO TO 102 
RETURN 

END 


A. 5 Subroutine FUNEV(K,X,Y,F) 


Description . This subroutine is called by RK7. It describes the 


system of differential equations for the different orders of the Taylor 


series. 


Parameters . K is the order of the Taylor series term. F is the 


matrix of the derivative values. Y is the matrix of the fundamental 


solution. W is the frequency. 


SUBROUTINE FUNEV (KiXiYiF) 

£ differential I eqn SUFPORTS ™ E RUNGE kutta and describes the 

c 

£ THI SSUBUILl- NEED CHANGING IF SEI IS TO HAVE DERIVATIVES 

£ IF THE TORSIONAL MODES ARE BEING DETERMINED 

£***♦ K=ORBER 0F THE TAYLOR SERIS TERMS 

C**»* F=DERIVATIVES VALUES 

C***» Y=FUND AMENT AL SOLUTION 

C***» N=FREQUENCY 

REAL*8 F» Y» DRM > U > SEI > X 
•DIMENSION F(8f13) »Y<8) 

COMMON/EGNVL/W 
F ( 1 . K ) =Y (2) 

F(2.K)=Y<3) 

F < 3 f K ) =Y < 4 ) 

F<4iK)=DRM(X)*l)**2*Y(l)/SEI(X) 

F < 5 f K ) =Y ( 6 ) 

F(6fK)=Y(7) 

F ( 7 , K ) =Y ( 8 ) 

F<8fK)=DRM < X)*U**2*Y(5) / SEI(X) 

RETURN 

END 
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A. 6 Function SEI(X) 


Description . This function determines the bending stiffness of the 
wing at any point of its span. SEI can be changed to fit any kind of 
wi ng . 

Parameters . X is the distance from the semi -span. 


FUNCTION SEI ( X ) 

C*** X=DISTANCE FROM SEMI SPAN 
IMPLICIT REAL*B (A-H.O-Z) 
SEI=900000000.DO 
IFIX.GE.ll .DO) SEI =90000000 .DO 
RETURN 
END 


A. 7 Su broutine SECANT(X1 ,F1 ,X2,F2,NC,NN) 

Description . This subroutine improves on the approximation given 
by BISECTION by decreasing the error bound on two successive values of 
the characteristic determinant. 

Parameters . XI and X2 are two consecutive values in the frequency 
domain; FI and F2 are the corresponding values of the characteristic 
determinant. EL and EU are the error bounds on the values of the calcu- 
lated solution of RK7. 


103 


500 


101 


SUBROUTINE SEC ANT < XI , F 1 . X2 , F2 , NC , NN > 

REAL*8 Xl»Fl»X2»F2fXhl»FP»DXiXM i FM » CCS i EU > EL 

EL =.00000000 IDO 

EU= . OOOOOOOIDO 

CCS= . 0001 DO 

XH1=X1 

CALL FUNtXl fFI » EL » EU > NC » NN ) 

CALL FUN(X2rF2>EL.EU.NC.NN) 

FP=<F2-F1 )/ CX2-X1 > 

DX=-F1/FP 

XM=X1+DX 

CALL FUN(XM»FM»EL»EU»NC»NN) 

IF< (DABS(XM-XH1 )/XM) .LT.CCS) GO TO 101 

IF (DABS(FM) .LT .CCS) GO TO 101 

IF(FM*F1. LE . 0 . DO ) GO TO 500 

X 1 =XM 

F 1 = FM 

XM1=XM 

GO TO 103 

X2=XM 

F2 = Fh 

Xhl =XM 

60 TO 103 

X2 = XM 

F2 = FM 

RETURN 

END 
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A. 8 Subroutine DMODE(YY,NN) 


Description . This subroutine calculates the values of the modes at 
151 points of the semi -span. These values are normalized such that the 
value of the mode at the wing tip is unity. 

Parameters . YY is the matrix of the values of the mode at NN = 151 
points of the semi -span. 


100 


SUBROUTINE DMODE<YY,NN) 

REAL*8 Y , X , Cl » C2 1 D 
COMMON/FACT/Y ( 8? 151 ) ,X( 151) 
DIMENSION YY < 151 ) 
»=Y<1'NN>*Y< 7 ,NN)-Y<5»NN>*Y<3»NN> 
Cl = Y < 7 , NN ) /D 
C2=-Y(3,NN)/D 
DO 100 1=1, NN 
YY<I)=C1*Y(1,I)+C2*Y(5,I) 

RETURN 

END 


A. 9 Function DRM(X) 

Description . This function gives the mass distribution along the 
semi -span as given in Figure 3. 

Parameters. X is the distance from the semi -span. 


FUNCTION DRM < X ) „ _ 
IMPLICIT REAL*8 <A-H,0-Z> 
DRM=130.D0 _ 

I F < X . LE . 4 . DO ) DRM=2205 . DO 
IF(X.LE.ll .DO. AND.X.GE.8.D0) 


RETURN 

END 


DRM=2600 . DO 


A. 10 Function RM(X) 

This function is the same as DRM(X). 


FUNCTION RM < X ) 

RM= 1 30 • 

IFCX.LE.4.) RM=2205. 

I F ( X . LE . 11. .AND.X.GE.8. ) RM = 2600. 

RETURN 

END 
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A.11 Function SIMP(Y.H.N) 


Description . This function is called by COF to integrate the 
functions to calculate the generalized masses and the aerodynamic cross 
terms of the modes. 

Parameters . Y is the mode matrix, H is the distance between two 
consecutive nodes, and N is the number of nodes. 


100 

200 


FUNCTION SIMP ( Y t H » N ) 

DIMENSION Y(N) 

T1 = 0 . 

Jl=N-2 
T2 = 0. 

DO 100 I =3 > J1 > 2 
T1=T1+Y< I ) 

DO 200 I =2 > N . 2 
T2=Y ( I ) +T2 

SIMP=H*< Y<1 >+Y<N)+2.*Tl+4.*T2)/6. 

RETURN 

END 


A. 12 Subroutine C0F(N,YY,M2,NN,I1 ) 

Description . This subroutine calculates the generalized masses and 
the aerodynamic cross terms of the modes. The rigid-body translation 
and rolling and the first four elastic modes are taken into account. 

Parameters . RP is the array of the six modes dealt with here. 
RP(1,I) represents the riqid-body translation and RP(6,I) represents 
the rigid-body rolling. A and B are matrices for the aerodynamic 
cross terms. RA is a matrix representing the semi -chord distribution 
along the semi-span. 


c 

c*** 

c**« 

c*** 

c*** 

c*** 

c*** 

c*** 

c*** 

C**» 


c*** 


601 


THIS^SUBROUT INE^SETS^UP^ THE 1 FUNCTIONS FOR INTEGRATION 
YY=NEW MODE THEN LATER USED AS SCRATCHED ARRAY 
RP= ARRAY OF MODES 
U=FREQUENCY 

GM = GENERALIZED MASS 

A X B=AREODYNAMMIC CROSS TERMS 

n=number P of N mode BEING UROKED ON 

NN=NUMBER OF NODES 
I 1=NUMBER OF TIMES COF CALLED 
C0HM0N/DATARA/RA<151) _ 

COMMON/TRAN 1/W( 6)jGM(6)»RP(6i151)tA(6>6>tB( 6»6) 
DIMENSION YY<151)»XX(151)iRR<151) 

READ DATA 
H=33,/(NN-1> 

IF(Il.GT.l) GO TO 500 
DO 601 I = 1 ? NN 
XX< I )=(I-1 )*H 
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600 

C*** 

500 

101 

102 


400 

103 

105 


108 


107 

104 


NM=N-1 

U<1)=0. 

GH ( 1 ) =40000 . 0 
DO 600 1=1 »NN 
RP< 1 » I)=l . 

A(l. 1>=45. 038925 
B(l,l)=52. 939763 
CALCULATE NEW DATA 
CONTINUE 
DO 101 1=1. NN 
RP ( N . I ) =YY < I ) 

DO 102 1=1. NN 
A1 =RH ( XX ( I > > 

B1=RP(N,I)*RP<N.I) 

YY< I )=A1*B1 
CONTINUE 

GH<N)=SIHP( YY»H»NN>*2. 

U(N)=W2 
DO 104 J=1»N 

IF(N.EQ.6)G0T0 105 
NR2=N/2. 

NR 1 = J/2 ♦ 

N2=N-2*NR2 

N1=J-2*NR1 

IF<N2.EQ.0.AND.N1 .EO.O) GO TO 400 

IF(J.EO.l) GO TO 400 

IF(N2.EQ.l.AND.Nl.EQ.l) GO TO 400 

A < N » J ) =0 . 0 

B<N. J)=0.0 

GO TO 105 

CONTINUE 

DO 103 1=1. NN 

YYd >=RA( I >*RP< J. I )*RP<N» I ) 

RR(I)=RA(I)*YY(I> 

A(N.J)=SIMP(RR.H.NN> 

B(N.J)=SIMP< YY.H.NN) 

CONTINUE 

I F ( N • NE . 6 ) GOTO 107 
A < 6 . 1 ) =0 . 0 
A (6. 2) =0.0 
A ( 6 . 4 ) =0 • 0 
IF(J.EQ.l) GOTO 107 
IF< J.EQ.2)G0T0 107 
IF ( J . EQ . 4 ) GOTO 107 
DO 108 1=1. NN 

YY(I)=RA<I)*RP< J»I)*RP<N.I> 
RR<I)=RA<I)*YY<I> 
A(6»J)=SIMP(RR.H.NN) 
B(6.J)=B1MP( YY.H.NN) 

CONTINUE 

CONTINUE 

RETURN 

END 


A. 13 Subroutine STRK(NN) 

Description . This subroutine sets up all the coefficients needed 
in RK7 to solve the system of differential equations using Runge-Kutta 
method of order seven. It also describes the semi -chord distribution 
along the semi -span. The semi -chord distribution is illustrated in 
Figure 4. 

Parameters . A, B, and C are the coefficient matrices. RA is a 
matrix representing the semi -chord distribution at 151 points. X is 
the abscissa matrix of the mode values. 
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SUBROUTINE STRK C NN > 

THIS SUBROUTINE INITIALIZES THE COEFICI 
FIXES THE NODES AND WING PLAN 
C*tt$ NN=NUHBER OF NODES 
C**»* AfBfC = COEFFICINTS FOR RK7 
Ct**t X=NODE VALUES 
Cttt* RA=UING PLAN ARRAY 

REAL*8 A,B,CfCHfXf Y,DEL 
COMMON/DATARA/RA( 151 ) 

CGHH0N/RKC/AC13) »B(13,12),C<13),CH<13) 
COHHON/FACT/YO, 151) ,X( 151 ) 

DEL =33 .DO/<NN-l) 

DEL2=33./(NN-1) 

ND=(NN-l)/3. 

RMl=<NN-l>/3. 

DO 100 I = 1 1 NN 
X(I)=(I-1)*DEL 
RA< I ) = 1 .0 

IF< I .GT.ND)RA<I)=1.-.027*<I-RM1)*DEL2 
100 CONTINUE 

A<1) =0 * ODO 

A(2)=2. DO/27. DO 

A < 3 ) = 1 . DO/ 9 .DO 

A(4)=1.D0/6.D0 

A(5)=5. DO/12. DO 

A ( 6 > = 1 .DO/2. DO 

A ( 7 > =5 . DO/6 . DO 

A<8)=1.D0/6.D0 

A < 9 ) =2 • DO/3 .DO 

A ( 1 0 ) = 1 * DO/3. DO 

A < 1 1 ) = 1 .DO 

A < 1 2 ) =0 . DO 

A<13)=1 .DO 

B ( 1 r 1 >=0. ODO 

B<2»1)=2. DO/27. DO 

B(3»l)=l. DO/36. DO 

B<3,2)=1 .DO/12. DO 

B< A f 1 > = 1 . DO/2 A . DO 

B ( A f 2 ) =0 . DO 

B(Af3)=1.D0/8.D0 

B < 5 f 1 ) =5 . DO/12 . DO 

B ( 5 , 2 ) = 0 * ODO 

B(5,3>=-25. DO/16. DO 

BIS, 4>=25. DO/16. DO 

B ( 6 » 1 ) = 1 • D0/20 . DO 

B ( 6 v 2 ) =0 . DO 

B(6»3)=0.D0 

B(6f A)= l .DO/A.DO 

B(6i5)=l .DO/5. DO 

BI7, 1>= -25. DO/ 108. DO 

B(7t2)=0.D0 

B( 7 f3 )=0 . ODO 

B(7. A)=125. D0/108. DO 

B(7»5)=-65. DO/27. DO 

B(7,6)=12S.DO/5A.DO 

B(8f1)=31. DO/300. DO 

B ( 8 f 2 ) =0 . ODO 

B(8f3)=0. DO 

B ( 8 f A ) =0 • DO 

B(8f5)=61 .DO/225. DO 

B(8f6)=-2.D0/9.D0 

B(8f7)=13. D0/900. DO 

B < 9 1 l ) =2 . DO 

B(9f2)=0.D0 

B ( 9 i 3 ) =0 . DO 

B<9»A)=-53. DO/6, DO 

B(9f5)=70A.D0/A5.D0 

B(9t6)=-107.DO/9.DO 

B(9.7)=67. D0/90. DO 

B < 9 f 8 ) =3 .DO 

B< 10 f 1 )=-91 . D0/108 .DO 

B(10»2)=0» ODO 

B(10f3)=0. DO 

B( 10, A) =23. DO/ 108. DO 

B< 10f5)=-976.D0/135.DO 

B( 10,6) =311 » BO/5 A • DO 

B<10,7>=— 19.DO/60.DO 

B(10,8)=17.D0/6. DO 

B(10,9)=-l .DO/12. DO 

Bill , 1) =2383. DO/ A 100. DO 

B< 1 1 , 2 ) =0 . DO 

B ( 1 1 , 3 ) =0 . DO 

B< 11 , A) = -3A1 .DO/ 16 A. DO 

B( 11,5 >=AA96. D0/1025. DO 

B(ll,6)=-301. DO/82. DO 

B(1 1,7) =2133. DO/A 100. DO 

B< 1 1,8)= A5. DO/82. DO 

B< 11 ,9) =A5.DO/16A,DO 

Bill, 10) =18. DO/41 .DO 

B(12,l)=3. DO/205. DO 

B<12,2>=0. DO 

B ( 12,3) =0 . DO 

B ( 12 , A ) =0 . DO 

B(12,5>=0. DO 
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B(12»6)= _ 6«D0/41. DO 
B< 12r7)=-3.D0/205.D0 
B<12>8)=-3. DO/41 .DO 
B(12»9>=3.DO/41.DO 
B(12»10)=6.D0/41.D0 
B( 12 » 1 1 ) =0 . ODO 
B ( 1 3 > 1 ) =-1777 • DO/41 00 . DO 
B ( 13 r 2 ) =0 . DO 
B ( 13 i 3 ) =0 • DO 
B( 13»4>=-341 .DO/ 164. DO 
B(13>5)=4496» DO/1025 . DO 
B113f6)=-2B9 . DO/82. DO 
B(13f7)=2193. DO/ 4 100 . DO 
B<13f8)=51 .DO/82. DO 
B<13»9> =33. DO/164. DO 
B(13f10>=12. DO/41. DO 

B ( 13 • 1 1 ) =0 . DO 
B( 13» 12 ) =1 • DO 
C(1 )=41 . DO/840 . DO 
C < 2 ) =0 , DO 
C < 3 ) =0 . DO 
C(4 ) =0 . DO 
C ( 5 ) =0 . DO 
C<6)=34 .DO/105. DO 
C<7>=9. DO/35. DO 
C(8)=9. DO/35. DO 
C(9)=9. D0/280. DO 
C( 10) =9. DO/280. DO 
C( 11 )=41 , D0/840 .DO 
C( 12 ) =0 . DO 
C ( 13) =0 . DO 
CH< 1 ) =0 . DO 
CH ( 2 ) =0 . DO 
CH(3) =0 .DO 
CH < 4 ) =0 . DO 
CHI 5 ) =0 . DO _ ^ 
CH(6)=34. D0/105. DO 
CH<7)=9. DO/35. DO 
CH(8)=9.D0/35.D0 
CH < 9 ) = 9 . B0/280 . DO 
CH(10)=9. DO/280. DO 
CHtll)=O.DO „ 

CH< 12) =41 .DO/840. DO 
CH ( 13 ) =41 • D0/840 . DO 
RETURN 
END 


A. 14 Subroutine RK7(NS,NN,EL,EU) 

Description . This subroutine is used to calculate the solution of 
the differential equation for a given frequency. It uses the Runge- 
Kutta-Fehlberg method seven and eight order. 

Parameters . EL and EU are lower bound and upper bound for the 
error. NS is the number system of the equations. Y is the solution 
matrix. 


c*** 

c**». 

c*** 

c*** 

c*** 

c*** 

c*** 


101 


SUBROUTINE RK7 < NS . NN . EL » EU ) 

RUNGE KUTTA FEHLBERG SEVENTH ORDER 
EL=ERROR LOUER BOUND 
EU=ERROR UPPER BOUND 
NS=NUMBER OF SYSTEM OF EON 
Y=SOLUTION 

NN-NUHBER OF PTS TO DETERMINE THE SOLUTION 
RL=LENGTH OF INTERVAL 
DIMENSION Y0<8) ,FI8. 13) , YY<8) 

DIMENSION DY4(8)#DY5(8)»Y1(8) 

do A ioi x^i Ins 5,xx ' YY,TH ' A,B,C,F,DD1 ’ cn2,H ’ Y1,YD ’ CH ’ EL ' EU,X ' y 

Y0<I)=Y(If1) 
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NT=NN-1 
L= 1 

C*** MAIN DO LOOP INCREAMENT TO EACH NODE 
DO 100 11=1 fNT 
NC=0 

H=X(I1+1)-X(I1) 

GO TO 203 
207 L=L-1 

GO TO 203 
204 L=L+1 

NC = 1 

203 DO 201 1=1 f NS 

201 Y 1 < I ) = YO ( I > 

TH=H/L 

C*** DO LOOP FOR STEPS BETWEEN NODES 
DO 200 12=1 »L 

C*** DETERMINE THE NEEDED FUNCTION EVALURATION 
DO 300 K=1 f 13 
KM=K-1 

DO 301 J= 1 v NS 

XX=X(I1)+TH*(I2-1)+ACK)*TH 
YY( J>=Y1< J) 

IF(KM.EQ.O) GO TO 303 
DO 302 13=1 r KM 

302 YY< J)=TH*B<Kf I3)*F< Jf I3)+YY< J> 

303 CONTINUE 
301 CONTINUE 

300 CALL FUNEV(KfXXfYYfF) 

C*** DETERMINE SOLUTION VALUE FOR END OF STEP 
DO 500 I = 1 f NS 
DY 4(1) =0 . 0 
500 DY5 ( I ) =0 . 0 

DO 401 1 = 1 f NS 
DO 402 K= 1 r 13 

DY4< I >=TH*C<K)*F< I fK)+DY4< I > 

402 DY5<I)=TH*CH(K)*F(IfK)+DY5<I) 

401 CONTINUE 

C*** ERROR AND STEP SIZE CONTROL 

DD1 = DABS< ( DY4 < 1 ) -DY5 ( 1 > ) / < DY4 < 1 > +Y 1 < 1 ) ) ) 
DD2=DABS( <DY4<5)-DY5<5> )/<DY4<5)+Yl <5) ) ) 
I F ( DD 1 . LT .EU.AND.DD2.LT . EU > GO TO 202 
GO TO 206 

202 IF (DD1 . GT . EL . AND . DD2 • 6T • EL ) GO TO 204 
IF(L.EQ.l) GO TO 204 

IF(NC.EQ.l) GO TO 204 
GO TO 207 

204 CONTINUE 

DO 205 1=1 f NS 

205 Y1(I)=Y1(I)+DY4(I) 

200 CONTINUE 

DO 102 1 = 1 f NS 
YO < I ) =Y1 < I ) 

102 Y<IfI1+1)=Y0(I) 

100 CONTINUE 
RETURN 
END 


A. 1 5 Subroutine D2 

Description . This program solves the forced vibration program. It 
calculates the amplitudes of the different modes due to sinusoidal gust 
at different stations of the wing and for different gust frequencies. 
First, it sets the input data and then it solves. the system of algebraic 
equations for different gust locations. 

Parameters . NN is the number of modes at which modes are calcu- 
lated and N2 is the number of gust locations. 
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SUBROUTINE B2 

C THIS PROGRAM DETERMINES THE AMPLITUDES OF THE DIFFERENT'" 

C HODES TO A SINUSODAL GUST AT THE DIFFERENT STATIONS ALONG 
C THE WING. SUB... DO TAKES CARE OF INPUT AND OUTPUT PLUS 

C SETS UP THE COEFFICIENTS THAT ARE DRIVING FREQUENCY INDEPENDENT. 

C SUBROUTINE COEF SET UP THE COEFFICIENT MATRIX FOR EACH DRIVING 

C FREQUENCY. WHILE SUB... GAUSS DOES HALF OF THE REDUCTION AND 

C SUB... BACKS FINISHES THE REDUCTION AND DOES BACK SUBSITUTION 
C FOR THE DIFFERENT NON-HOMOGENOUS VECTORS CORRESPONDING TO 
C DIFFERENT GUST LOCATIONS. 

C*** NN=NUMBER OF NODES THE MODES ARE DETERMINE ON 
C*** N2=NUMBER OF GUST LOCATION 
NN=151 
N2 = 20 

CALL DO ( NN » N2 ) 

RETURN 

END 


A. 16 Subroutine D0(NN,N2) 

Description . This subroutine sets the constants needed later in 
calculations of the mode amplitudes. Then it normalizes the generalized 
masses and the aerodynamic cross terms, after which different gust 
frequencies are considered. For each of these frequencies, the coeffi- 
cient matrix of the unknowns and the nonhomogeneous vectors for 19 gust 
stations along the semi-span are calculated by calling subroutine COEF. 
Subroutine GAUSS is called to perform reduction of the coefficient's 
matrix of the unknowns. Finally, subroutine BACKS is called to perform 
reduction on the nonhomogeneous vectors and to perform back substitution. 

Parameters . Y is the solution of the system of equations, i.e., 
the amplitudes of modes. RP is the array of the six modes considered. 

X is the gust location. A and B were arrays of the aerodynamic cross 
terms and now become the arrays of the normalized aerodynamic cross 
terms. W is the array of the natural frequencies of modes. GM was the 
array of the generalized masses of the modes; these values are normal- 
ized and transferred to the array GAMA. OMEG is the reduced natural 
frequency; these frequencies are reduced with respect to the semi -root 
chord BR and flight speed of the airplane U = 575 ft/sec. RRK(37) is 
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the matrix of reduced frequencies of the gust. RO = 0.0765 is the air 

2 

density. S = 960 ft is the wing area. 


c 

c 

c 

c*** 

c*** 

c*** 

c*** 

c*** 

c*** 

c*** 

c**t 


c *** 


602 

601 


100 


2 

1 

102 

500 


SUBROUTINE D0<NN.N2> 

THIS SUBROUTINE TAKES CARE OF INPUT ANB OUTPUT 
AND PERFORMS OPERATIONS THAT ARE DRIVING FREQUENCY 
INDEPENDENT, MEANING GAMA AND OMEG. 

Y=SOLUT I ON » AMPLITUDES OF MODES 
RP=MODE ARRAY 
X=GUST LOCATION 

A, B=AERO DYNAMIC CROSS PRODUCTS 
U=NATURAL FREQUENCIES OF MODES 
GM=GENERALIZED MASS OF HODES 
OMEG=REDUCED NATURAL FREQUENCY 
GAMA=NON DIMENSIONAL GM 
COMMON/SOL/Y (20,6) 

COMMON/TRAN 1/U( 6) ,GM< 6) ,RP <6,151 1 ,A<6,6>,B<6,6> 
C0MH0N/TRAN2/RRK(37> 

C0MM0N/TRAN3/SY< 37 ,19,6) 

COMMON/EDAT/EEE 
COMMON/FAC/F I » BRiS.RO.U 
DIMENSION X ( 20 > 

COMPLEX Y,SY 
COMPLEX CMPLX 
C0MM0N/DAT/GAMA(6) , OMEG ( 6 ) 

N = 37 

PI=3. 14159 
BR=19. 0/2.0 
S=960.0 
R0= • 0765 
U=575.0 

PERFORM ARITHMATIC 
DO 601 1=1,6 

GAMAd >=GM(I )/<PI*RO*S*BR) 

OMEGC I )=W< I ) *BR/U 
JN = I 

DO 602 J=1»JN 
A(I , J)=A( I , J)*BR/S 
Bd, J)=B(I, J)*BR/S 
A(J,I)=Ad,J) 

B < J , I ) =B < 1 1 J ) 

CONTINUE 
CONTINUE 
DO 500 1=1,37 
IFd.LE.10) RK = I/1000. 

IFd.LE.lV.AND.I.GT.10) RK= ( 1-9 )/100. 

IFd.LE. 28. AND. I .GT.19) RK= (I -181/10. 

IFd.GT.28) RK=(I-27) 

DEL=33./(N2-1 > 

N2M=N2-1 

DO 100 11=1, N2M 

Xdl)=DEL*dl-.5> 

CALL C0EF(RK,X,DEL,N2M,NN) 

CALL GAUSS ( 6 , N2M ) 

CALL BACKS < N2M , 6 > 

RRK (1 1 =RK 

DO 1 J= 1 , N2M 

DO 2 J2= 1 , 6 

SY< I , J, J2)=Y( J, J2) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

RETURN 

END 


A. 17 Subroutine C0EF(RK,X,DEL,N2M,NN) 

Description . This subroutine calculates the matrix of coefficients 

of the system of linear algebraic equations. The coefficients of the 

homogeneous system are calculated only once since they do not change for 

all the gust locations. The nonhomogeneous vector has to be calculated 

for 38 gust stations along the entire wing. 
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Parameters . C is the coefficient matrix of the homogeneous system. 
D is the array of the nonhomogeneous vectors. RK is the reduced fre- 
quency of the gust. X is the matrix of locations of the gust. DEL is 
the distance between two gust locations. N2M is the number of gust 
locations. 


SUBROUTINE COEF ( RK » X . BEL , N2M , NN ) 

C THIS SUBROUTINE SETS UR COEFFICIENT MATRIX AND THE DIFFERENT 
C NON HOMOGENOUS VECTORS. 

C** C=COEFFICIENT MATRIX 

C*¥ D=ARR A Y OF NON HOMOGENOUS VECTOR 

C*** C=COEF IC I ENT ARRAY 

c*** RK=REDUCED FREQUENCY 

C**¥ X=LOCAT I ON OF OUST 

CGMM0N/LS/C(6»6) .IK 20 .6) 

COMMON/ DAT/GAMA (6) . OMEG ( 6 ) 

COMMON/TRAN 1/U(6).GM(6) . RP < 6 . 1 5 1 ) . A < 6 . 6 ) r B ( 6 . 6 > 
COMMON/FAC/PI .BR.S.RO.U 
DIMENSION X C 20) . SI <6 . 6) 

COMPLEX C.CI.D.CC.RKK 
COMPLEX CMPLX 
C*** READ DATA 

CI=CMPLX(0. . 1 . > 

DO 103 1=1.6 
DO 104 J=l,6 

C(I» J>=-RK¥*2*A(I, J)+2*CI*RK*CC(RK>¥B(I. J) 

IF( I . NE . J) GO TO 901 

C(I. J>=C(I . J)+GAMA(I)*(OMEG(I )**2-RK**2> 
S1<I.J)=CABS(C(I,J>) 

CONTINUE 
CONTINUE 
GOTO 13 
CONTINUE 
DO 106 J= 1 . N2M 
DO 105 1=1.6 

D< J. I)=2*BR/S*RA<XCJ) ) tRPH (I.X(J). NN ) ¥RKIK C RK ) ¥DEL 
CONTINUE 
CONTINUE 
RETURN 
END 


901 

104 

103 

13 


105 

106 


A. 18 Complex Function CC(RK) 

Description . This function calculates the Theodorsen function in 
terms of the Bessel functions RJO, RJ1 , RYO, and RY1 . 

Parameters. CC(RK) is a function of reduced gust frequency RK. 


COMPLEX FUNCTION CC<RK> 

C THIS FUNCTION CALCULATES THE THEODORSEN FUNCTIONS 
COMPLEX Cl 
COMPLEX CMPLX 
CI=CMPLX(0. .1 . ) 

PJ1 =R J1 (RK) 

P JO=R JO ( RK ) 

PY1=RY1 (RK) 

P YO=R YO ( RK ) 

F=FJ1¥(PJ1+PYO)+PY1*(PY1-PYO) 

G=PY1*PY0+PJ1*PJ0 

H=(PJ1+PY0)**2+(PY1-PJ0)**2 

CC=(F+CI*G)/H 

RETURN 

END 
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A. 19 Function RPH(I,Y,N) 

Description . This function calculates an average value for the 
modes stored in RP at a certain gust location. 

Parameters . I is the mode number, Y is the gust location, and N is 
the number of modes. 


500 


600 


FUNCTION RPH ( I , Y > N ) 

COMMON/TRAN1/W(6) >GH(6) ,RP(6»151> ,A(6,6> ,B(6f6) 
DEL— 33*/(N— 1) 

NN=ABS<Y)/DEL+1 
RS=ABS < Y > /DEL+1 . -NN 

RPH=RP(IfNN)+RS*<RP<IfNN+1)-RP<IfNN) )/DEL 
IF(Y.LT.O.O) GO TO 500 
GO TO 600 
CONTINUE 

IF<I .EQ.3.0R. I .EQ.5) RPH=-RPH 
IF(I.EQ.6)RPH=-RPH 
RETURN 
END 


A. 20 Function RA(Y) 

Description . This function describes the semi -chord distribution 
along the semi -span. 

Parameters . Y is the distance from mid-span. 


FUNCTION RA ( Y > 
RA=1-.027*<ABS<Y>-11. ) 
IF(ABS(Y) .LE.ll) RA~1 • 
RETURN 
END 


A. 21 Function RY1 (X) 

Description . This function is the Bessel function of the second 
kind order one. 

Parameters. X is the reduced frequency of the gust. 


FUNCTION RY1(X) 

Z= ( X/3 . >**2 

RY1 = < < < ( ( < . 0027873*Z-. 0400976 >*Z+. 3123951 ) *2-1 .31 64827 )*Z 
++2 . 1682709 ) *Z+ . 221209 1 ) *Z- .6366198+ . 6366198*X*AL0G( X/2 . ) *R J1 ( X ) 
+ )/X 
RETURN 
END 
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A. 22 Complex Function RKK(RK) 


Description . This function is the Kussner function. It is calcu 
lated from Theodorsen function CC(RK) and Bessel functions RJO and RJ1 
Parameters . RKK(RK) is a function of reduced gust frequency RK. 


COMPLEX FUNCTION RKMRK) 

C FUNCTION DETERMINES THE GUST FORCE FUNCTION 
COMPLEX CI.CC 
COMPLEX CMPLX 
Cl = CMPLX ( 0 . >1. ) 

PJ 1=R J 1 ( RK > 

RKK=CC(RK)*<RJO<RK)-CI*PJl )+CI*PJl 

RETURN 

END 


A. 23 Function RJ1 (X) 

Description . This function is the Bessel function of first kind 
order one. 

Parameters. RJ1 is a function of reduced gust frequency RK. 


FUNCTION RJl(X) 

Z=<X/3. >**2 

RJ1=<<<<<<. OOOO 1109CZ-. 00031761) *Z+ . 00443319 ) *Z- . 039542B9 > *Z 
++. 2 1093573 >*Z~. 56249985 )*Z+ .5>*X 
RETURN 
END 


A. 24 Function RJO(X) 

Description . This function is the Bessel function of first kind 
order zero. 

Parameters . RJO is a function of reduced gust frequency RK. 


FUNCTION RJO ( X ) 

+Z- 2 T 2499997 )*Z+1 Z 0* 0039444 ) * z+ ‘ 0444479 ) *Z-. 31 63866 ) *Z+1 . 2656208 >* 

RETURN 

END 


A. 25 Function RYO(X) 

Description . This function is the Bessel function of second kind 
of order zero. 


Parameters . RYO is a function of reduced gust frequency RK. 
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FUNCTION RYO ( X > 

Z=<X/3. >**2 

RYO=( ( < ( ( -.00024846*Z+. 0042791 6 )*Z- .04261214 >*Z+. 253001 17) *Z 

+ orT.125 0384) * Z+ * 60936 > * Z+ ’ 3A746691 +.63A619B*RJ0CX ) *AL0G(X/2. ) 
RETURN 
END 


A. 26 Subroutine GAUSS(N,N2) 

Description . This subroutine performs Gaussian elimination on the 
coefficients matrix of the homogeneous system. 

Parameters . N is the number of equations and N2 is the number of 


c*** 

c*** 

c*** 

c*** 


104 

103 


105 


102 

101 

100 


SUBROUTINE GAUSS(NfN2) 

THIS SUBROUTINE DOES GAUSSIAN ELIMINATION FOR ONLY THE COEFFICIENT 
MATRIX. SCALED PARTIAL PIVOTING IS USED. 

COMM0N/LS/C<6fA)fD(20fA) 

COMMON/P IVOT/ I PEN< A) 

DIMENSION S < 6 ) 

COMPLEX CfD 
K=P I VOT INDEX 
C = COFF I C I EN T ARRAY 
D=INHOMOGENOUS VECTOR 
N=NUMBER OF EQN 
DO 103 I = 1 » N 
IPEN( I )=I 
S<I)=0. 

DO 104 J= 1 f N 

IF(CABS(C<IfJ)).GT.S(I>) S<I)=CABS<C<IfJ)) 

CONTINUE 

NM=N-1 

DO 100 KK=1fNM 
I S=KK+ 1 
IF=IFEN(KK> 
j — K k 

CM=CABS<C(IPfKK) )/S(IP) 

DO 105 I=ISfN 
I P= I PEN ( I > 

T=CABS(C< IP.KK) )/S( IP) 

IF(T.LE.CM) GO TO 105 

CM = T 

J=I 

CONTINUE 
I F'K = I PEN ( J ) 

IF'EN< J) = IPEN(KK) 

I F'EN < KK ) =IPK 
DO 101 II=ISfN 
I = I PEN (II) 

K=IPEN<KK) 

C ( I f KK ) =C ( I f KK ) /C ( K f KK ) 

DO 102 J=ISfN 

C( I f J)=C(If J)-C< I fKK)*C(Kf J) 

CONTINUE 

CONTINUE 

RETURN 

END 


A. 27 Subroutine BACKS(N1,N2) 

Description . This subroutine performs reduction on all the non- 
homogeneous vectors corresponding to the gust locations on the semi- 
span. Then it makes back substitution and calculates the unknowns. 
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Parameters . NT = 19 is the number of the nonhomogeneous vectors 
each corresponding to a gust location along the semi -span. N2 = 6 is 
the dimension of the nonhomogeneous vectors. 


SUBROUTINE BACKS<N1,N2> 

C DOES REDUCTION ON THE NON HOMOGENOUS VECTOR AND THEN DOES 
C BACK SUBSI TUT I ON . 

C*** N 1 =NUMBER OF NON HOMOGENOUS VECTOR 
C*** N2=D I MENSI ON OF NON HOMOGENOUS VECTOR 
C0HM0N/PIV0T/IFEN<6) 

COMM ON /L S/C <6>6) i D < 20 > 6 ) 

COMMON/SOL/Y < 20 f 6 > 

COMPLEX C f D f Y 
C*tt Kl— SOLUTION INDEX 

C*** REDUCTION ON NON HOMOGENOUS VECTOR 
DO 100 Kl-1 f N1 
I P = I PEN < 1 ) 

Y<K1 f 1 )=D<K1 , IP) 

DO 101 KK=2fN2 
K = IF'EN<KK) 

T = 0.0 
JN=KK- 1 
DO 102 J= 1 v JN 

102 T=C(K, J)*Y(K1 f J>+T 
101 Y<K1 fKK)=D<K1»K)-T 

Y(K1fN2)=Y(K1 fN2)/C<KfN2) 

C*** BACK SUBSI TUTI ON 
JJ=N2 

DO 103 K=2»N2 
JS = JJ 
JJ=JJ-1 
KK=IPEN( JJ) 

T = 0 . 0 

DO 104 J=JS.N2 
104 T=C(KKf J)*Y(K1» J)+T 

103 Y<K1fJJ)=<Y<K1fJJ)-T)/C(KKfJJ) 

100 CONTINUE 

RETURN 

END 


A. 28 Subroutine D3 


Description . This program calculates the power spectrum of the 


wing tip velocity. After assigning a turbulence scale it calls subrou- 


tine SPEC to calculate the output spectrum corresponding to this scale. 


Parameters . N2 is the number of gust locations and N is the number 


of gust frequencies. TL is the turbulence length scale. 


SUBROUTINE D3 

C THIS PROGRAM PERFORMES THE ARTHIMETIC TO DETERMINE WING TIP 

C VELOCITY POWER SPECTRUM. SUB... SPEC PERFORMS THE CALCULATION 

C AND FUN ... TSPEC EVALUATES ATMOSPHERIC TURBULENCE SPECTRUM. 

C*** N2 = NU(1BER OF GUST STATIONS 
C*** N=NUMBER FOR DRIVING FREQUENCIES 
C*** TL= TURBULENCE LENGTH SCALE 
COMMON/T RAN2/RRK < 37 ) 

COMMON /TRAN3/SY <37 * 19»6) 

COMMON/E D AT /EEE 
COMPLEX SY 
CALL C0EF1 
CALL C0EF2 
BR=19./2. 

N=37 

N22=20 

N22=N22-1 
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N2=N22*2 
TL=21 12 . 

DO 100 1=1 , N 
RK=RRK< I ) 

RNU=RK*TL/BR 
SS = 0.0 
U=575 . 

TS=TSPEC<SS,RNU,U,TL) 
IF ( TL . EQ • 132. )GOTO 

1000 


IF ( TL * EQ • 660 . IGOTO 

1001 


IF (TL.EQ.2112. JGOTO 

1000 

1000 

WRITE<30»2020) RK , TS 


1001 

GOTO 3060 

WRITE < 31 ,2020) RK.TS 


1002 

GOTO 3060 

URITE(40,2020)RK,TS 


3060 

CONTINUE 


2020 

FORMAT (2E13.6) 



CALL SPEC ( RK,RR,TL,N2 

, I , RP > 

100 

CONTINUE 


101 

CONTINUE 



RETURN 

END 



A. 29 Subroutine SPEC(RK,RR,TL,N2,IC,RP) 

Description . This subroutine determines the output spectrum. The 
cross spectrum between the output and the input is also calculated. 

Parameters . RR is the output spectrum; RK is the reduced natural 
frequency, RK = WxBR/V; Z is the response to the gust at one station; 

SS is the separation distance nondimensionalized by the turbulence scale 
RNU is a reduced gust frequency, RNU = RK*TL/BR where RK is the gust 
frequency. 


c 

c 

c*** 

c*** 

c*** 

c**» 


c»*» 


100 

c*** 


101 


101 

c*** 


SUBROUTINE SPEC < RK , RR , TL . N2 , IC . RP > _ 

THIS SUBROUTINE DETERMINES THE SPECTRUM OF THE WING 
TIP VELOCITY. 

RR= TOTAL AIRPLANE RESPONSE 
RK=REDUCED FREQUENCY UB/U 
2= RESPONSE TO GUST AT ONE STATION 
N2=NUHFER OF GUST STATIONS 
COMMON /TRAN 2 /RRK ( 37 > 

C0HN0N/TRAN3/SY ( 37 * 19 , 6 ) 

DIMENSION YC20.6) ,Z<40> 

COMPLEX Y,CI»Z»T»T2,SY»H»Z1»Z2»H1,H2,Z3»Z4 
READ DATA 
C1=CMPLX<0. >1.) 

BR=19./2. 

N22=N2/2 

DEL=66./N2 

U=575. 

DO 100 J = 1 , N22 

DO 2 J2=l r 6 

Y(J» J2)=SY(IC»J.J2) 

CONTINUE 

CONTINUE 

DETERMINE PLANES RESPONSE 
DO 101 J=1»N22 
T = 0 . 0 
T2=0 . 0 

DO 102 1=2.5 
T2 = Y < J > I > +T2 

IF< I . EQ.3.0R.I ,E0.5» T2=T2-2 . *Y ( J . I ) 

T=Y(J.I)+T 
Z( J)=T2*RK*CI 
Z( J+N22)=CI*RK*T 
DETERMINE PLANES TOTAL RESONSE 
TT=0 . 0 

DO 300 1=1, N2 
IS=I-1 
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T=0.0 

JN=N2-IS 

DO 301 J=1.JN 

301 T=2*REAL(Z< J)*C0NJG<Z< J+IS) ) >+T 
SS=IS*DEL/TL 
RNU=RK*TL/BR 
IF(IS.E0.0)T=T/2. 

300 TT=TSPEC(SS.RNU.U.TL)*T+TT 
RR=TT 

IF(TL.EQ. 132. )G0TD 7 
IF < TL . EQ . 640 . >G0T0 8 
I F < TL . EQ . 21 1 2 . > GOTO 7 

7 URITE(31»17644)RK.RR 
GOTO 3040 

8 URITE<20. 176441RK.RR 
GOTO 3040 

9 WRITE<41»17644) RK » RR 
3040 CONTINUE 

C*** DETERMINE PHASE OF PLANES RESPONSE 
T = 0 . 0 

DO 600 J=1»N2 
IS=N2-J 
SS=IS*DEL/TL 
RNU=RK*TL/BR 

T=TSFEC(SS.RNU»U.TL)*2*REAL(Z< J> >+T 


600 

CONTINUE 
RP = T 

IF<TL.EQ,132*)GOTO 17 
IF(TL.EQ.660.)G0TO 18 
IF(TL.EG.2112. )GOTO 17 

17644 

FORMAT <2E13»6) 

17 

WRITE (32. 17644) RK . RP 
GOTO 3050 

18 

URITE<7. 17644) RK . T 
GOTO 3050 

19 

WRITEC 42? 17644)RK»T 

3050 

CONTINUE 


RETURN 

END 


A. 30 Function TSPEC(SS ,RNU ,U,TL) 


Description . This function calculates power spectrum and cross- 


power spectrum from von Karman spectrum functions. This function is 


called in SPEC in order to calculate the output and the cross spectrum 


between the input and the output. 


Parameters . SS is a separation distance divided by turbulence 


scale. RNU = WxTL/BR is the reduced frequency of turbulence. TL is 


the turbulence length scale. 


FUNCTION TSPEC(SS.RNU»U.TL> 

C THIS FUNCTION DETERMINES TURBULENCE CROSS AND POWER 
C SPECTRUM FROM THE VON KARMAN SPECTRUM FUNCTION. 

C*** SS = SEF*ERATION DIVIDED BY TL < TURBLUNCE LENGTH SCALE) 

C*** RNU = W*TL/U THE REDUCED FRQUENC Y OF TURBULENCE 
C*** U=FL I GHT SPEED OR MEAN WIND SPEED 
C*** TL=TURBULENCE LENGTH SCALE 
IF(SS.EO.O.O) GO TO 500 
C CROSS SPECTRUM 

Z=SS*SQRT< 1 . +< 1 . 339*RNU)**2)/1 .339 

TSPEC = TL*. 10653/U*(4. 781 12*SS*» <5 ./3. >/Z**C5./6. )*BSL1 (Z) 
+-SS**( 11 ./3.)/Z**<ll./6.)*BSL2<Z») 

RETURN 

500 CONTINUE 
C POWER SPECTRUM 

Z=< 1 .339*RNU>**2 

TSPEC=TL*( l+<8./3. )*Z)/( 1+Z)**( 11 ,/6 . )/3. 14159/U 

RETURN 

END 
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A. 31 


Subroutine C0EF1 


Description . This subroutine is needed in order to calculate the 
coefficients of the polynomial that approximates the modified Bessel 
function of the second kind of order 5/6. 


c 

c 

c 


100 


101 


200 


SUBROUTINE C0EF1 

THIS SUBROUTINE SETS UP THE COEFF I ENETS FOR THE POLYNOMIAL 
APPROXIMATION FOR THE MOBIFRIEIi BESSEL FUNCTION OF THE 
SECOND KIND 5/A ORDER. 

C0HM0N/K13/A< 10) . B ( 10 > » A2 ( 10 » 

F=5./6. 

A( 1 ) = 1 .0/ .9405612296 
DO 100 1 = 1,9 
A(I+1)=A(I)/I/CF+I> 

F=1 . 0-F 

B(1 )=1 .0/5.56756615 

DO 101 1 = 1,9 

B< 1+1 >=B(I )/I/(F+I— 1 .0) 

S=4.*<5./6. )**2 
A2< 1 ) = 1 . 


DO 200 I = 1 t 9 

A2(I+l)=A2<I>*(S-<2*I-l>**2)/8./I 

RETURN 

END 


A. 32 Function BSL1(Z) 


Description . This function calculates the modified Bessel function 


of the second kind of order 5/6. 

Parameters . Z is a function of the reduced frequency of the turbu- 
lence, z = (SS/1 .339)/! + (1 . 339RNU) 2 . 


FUNCTION BSLl(Z) 

C THIS FUNCTION EVALUATES THE MODIFRIED BESSEL FUNCTION 
C OF THE SECOND KIND 5/6 ORDER 

C0MM0N/K13/A( 10)iB<10)«A2(10) 

IF < Z • LE . 2 ) GO TO 100 
Y=1 . /Z 

BSL1 = SCIRT< 1 .5707*Y)*EXP(-Z)*P0LY< A2, 10. Y> 

RETURN 

100 Y=(Z/2.0>**2.0 

RIP=(Z/2.0>**<5./6. >*POLY(A, 10, Y) 

RIN=POL Y(B. 10. Y>/< ( Z/2 . 0 ) ** ( 5 . /6 . > ) 

BSL1= <3. 141/2/SIN <5.0*3. 141/6.0) )*<RIN-RIP> 

RETURN 

END 


A. 33 Function POLY(A,N,Z) 


Description . This function performs polynomial evaluations. These 
polynomials are used in BSL1 or BSL2. The coefficients of these poly- 
nomials are calculated in C0EF1 and C0EF2. 
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Parameters . A is the matrix of coefficients of the polynomial. N 

o 

is the number of terms in the polynomial. Y is either 1/2 or (Z/2) 
depending on whether Z is less than or greater than 2. 


FUNCTION POLY<A.N,Z> 

: THIS FUNCTION DOES THE POLYNOMIAL EVALUATIONS 

DIMENSION AtN > 

T=A(N)*Z 

NN=N-2 

DO 100 1 = 1, NN 
100 T=< T+AIN-I ) >*Z 

P0LY = T + A(1 ) 

RETURN 

END 


A. 34 Function BSL2(Z) 

Description . This function evaluates the modified Bessel function 
of the second kind of order 11/6. 

Parameters . Z is a function of the reduced frequency of the turbu- 
lence, Z = (SS/1.339)/! + (1 .339RNU) 2 . 


THIS N FuJ°TION L EVALUATES THE MODIFIED BESSEL FUNCTION OF THE 
SECOND KIND 11/6 ORDER. 

COMMON/K23/E <10),G<10),E2<10> 

IF ( Z . LE . 2 ) GO TO 100 

BSL2=SQRT<1.5707*Y>*EXP<-Z>*P0LY<E2,10,Y> 

RETURN 

° RlP^f Z/?! 0*** < 1 1 . /6. ) *POLY<E, 10, Y) 

RIN = P0LY(G,10.Y>/< <Z/2.0)**<ll./6. > 

BSL' I =(3. 141/2/ SI N< 11. 0*3. 141/6. 0))*(RIN-R1F> 

RETURN 

END 


A. 35 Subroutine C0EF2 

Description . This subroutine calculates the coefficients for the 
polynomial approximation of the modified Bessel function of the second 
kind of order 11/6. 


SUBROUTINE C0EF2 

C THIS SUBROUTINE SETS UP THE COEFFICIENTS FOR THE POLYNOMIAL 
C APPROXIMATIONS OF THE MODIFIED BESSEL FUNCTION OF THE 
C SECOND KIND 11/6 ORDER. 

COMHON/K23/EUO) ,G< 10) »E2<10) 

c***** F ~one / over the gamma value of 1+ORDER ******* 

E<1>=1. 0/1. 724362254 
DO 100 1=1,9 

100 E(I+1>=E<I>/I/(F+I> , 

C***** ONE MINUS THE ORDER OF THE MODFIED BESSEL ** 

F= 1 • 0-F 
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101 

200 

* 


§o 1 io}-?^7f- 68107938> 

E2<1)=1. 

DO 200 1=1,9 

rIturn >=E2(I, * <s " (2 * i ~ 1) ** 2)/b * /i 

END 
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